Binary fluids under steady shear in three dimensions 
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We simulate by lattice Boltzmann the steady shearing of a binary fluid mixture with full hydro- 
dynamics in three dimensions. Contrary to some theoretical scenarios, a dynamical steady state 
is attained with finite correlation lengths in all three spatial directions. Using large simulations 
we obtain at moderately high Reynolds numbers apparent scaling exponents comparable to those 
found by us previously in 2D. However, in 3D there may be a crossover to different behavior at low 
Reynolds number: accessing this regime requires even larger computational resource than used here. 

PACS numbers: 64.75.-|-g, 47.11. Qr 
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Systems that are not in thermal equilibrium play a cen- 
tral role in modern statistical physics [Q. They include 
two important classes: those evolving towards Boltz- 
mann equilibrium (e.g., by phase separation following a 
temperature quench), and those maintained in nonequi- 
librium by continuous driving (such as a shear flow). Of 
fundamental interest, and surprising physical subtlety, 
are systems combining both features — such as a bi- 
nary fluid undergoing phase separation in the presence 
of shear. Here a central issue [|[ |^ is whether coars- 
ening continues indefinitely, as it does without shear, 
or whether a noncquilibrium steady state (NESS) is 
reached, in which the characteristic length scales Lx,y,z 
of the fluid domain structure attain finite 7-dependent 
values at late times. (We define the mean velocity as 
Ux = iy so that a;, y, z are velocity, velocity gradient and 
vorticity directions respectively; 7 is the shear rate.) 

Our recent simulations, building on earlier work of oth- 
ers ||], have shown that in two dimensions (2D), a 
NESS is indeed achieved [||. In 3D, the situation is 
more subtle. Fourier components of the composition field 
whose wavevectors lie along the vorticity direction feel no 
direct effect of the mean advective velocity [|[ [tJ . There- 
fore it might be possible for coarsening to proceed indef- 
initely by pumping through tubes of fluid oriented along 
z Q. Another crucial difference is that in 2D fluid bi- 
continuity is possible only by fine tuning to a percolation 
threshold at 50:50 composition (assuming fluids of equal 
viscosity) so that the generic situation is one of droplets. 
(Indeed, for topological reasons, droplets are implicated 
even at threshold Q.) In contrast, in 3D both fluids re- 
main continuously connected across the sample through- 
out a broad composition window either side of 50:50. 

In 3D experiments, saturating length scales are report- 
edly reached after a period of anisotropic domain growth 
However, the extreme elongation of domains along 
the flow direction means that, even in experiments, finite 
size effects could play a role in such saturation Q . Theo- 
ries in which the velocity does not fiuctuate, but does ad- 
vect the diffusive fiuctuations of the concentration field, 
predict instead indefinite coarsening, with length scales 
Ly^z scaling as 7-independent powers of the time t since 



quench, and (typically) Lx ^ jtLy ||9[]. As emphasized 
in 1^, in real fiuids, however, the velocity fiuctuates non- 
linearly in response to the advected concentration field, 
and hydrodynamic scaling arguments, balancing interfa- 
cial and either viscous or inertial effects, predict satura- 
tion instead e.g., L/Lq ~ (7To)-i or L/Lq ~ {jTo)'^/^ 
|, |ll|. Here, Lq = i^^/ipcr) and Tq = iy'^/{pa^), with 
p density, v = rj/p kinematic viscosity and a interfacial 
tension, are the characteristic length and time at which 
inertial effects start to influence coarsening Given 
these uncertainties as to the fate of sheared binary fluids 
in 3D, computer simulations of such systems, with full 
hydrodynamic velocity fluctuations, arc of great interest. 

Such simulations also offer demanding challenges to the 
state of the art in computational physics. The 2D lat- 
tice Boltzmann (LB) results of Q were obtained from 16 
production runs involving lattices ranging from 512 x 256 
to 2048 X 1024 (all systems having aspect ratio 2:1). 
Many pre-production runs were required to steer sim- 
ulation parameters so as to avoid finite-size effects and 
other artefacts. This effort was rewarded, however: the 
unique parametric fiexibility of LB allowed us to probe 
over six decades of reduced shear rate jTo Q. Below, 
we extend that work to three dimensions with 9 produc- 
tion runs on 512 x 256 x 256 lattices, and 3 larger runs 
of 1024 x 512 X 512 (i.e., aU with aspect ratio 2:1:1). 
Even given the excellent parallel scaling of LB on multi- 
processor machines, each one of these 12 datasets re- 
quired more computational resource than the entirety of 
Ref. 1^. The production runs reported here were per- 
formed using 1024 processors of the IBM Blue Gene/L 
machine at the University of Edinburgh. 

Although our simulations arc not the first to address 
sheared binary fluids in 3D (see e.g. |l3)), earlier stud- 
ies have offered only inconclusive evidence of NESS for- 
mation in systems free of finite size effects. Such effects 
can cause fully lamellar or hexagonal cylindrical domains, 
which wrap the periodic boundary conditions with simple 
topologies that prevent further hydrodynamic coarsening 
14 ; but this "trivial" route to NESS relies directly on 
the periodic boundary conditions and is thus not avail- 
able in the bulk-system limit. Below we present evi- 
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Figure 1: Snapshots of the interface position at 7T0 — 22.47 
(top) and 47.45 (bottom) witli parameter set R019 (Table 1). 
These are representative of the observed NESS. The mean 
flow is rightward along the upper face of the simulation box 
and leftward at the lower face; the line of sight lies close to the 
vorticity (neutral) direction, z. Yellow and blue isosurfaces 
are constructed at <j) — ±0.2 to create a dividing surface color- 
coded by the adjacent fluids (both shown transparent). 



dence of NESS formation in systems retaining the com- 
plex topology expected in bulk samples, where a steady 
state dynamical balance can arise between the coarsening 
of bicontinuous domains under the action of interfacial 
tension, and their stretching by the flow (Figj^). 

The required parameter steering would not have 
proven possible without having the 2D runs to initially 
guide our selection — a methodology that can only suc- 
ceed if the physics in 2D and 3D is not radically differ- 
ent. Below we find that to be true for the upper few 
decades of the range of (^yTo)"^ addressed in within 
this range, evidence is given below for saturation of cor- 
relation lengths with finite values in all three directions. 
We then combine datasets using a quantitative scaling 
methodology developed for the unsheared problem in |12f| 
and for shear in |^; this allows scaling exponents to be 
estimated using combined multi-decade fits. Caution is 
required here due to residual finite size effects; these are 
unavoidable, particularly at high shear rates where we 
find NESS hardest to achieve numerically. Note that 
high shear rates correspond to low Reynolds numbers Re 
~ Ly'y/h' (due to the decrease of domain size with shear 
rate) ; these results could therefore signify new physics at 
low Rc B. However, much larger systems sizes might be 



needed to gain full access to this regime. 

The governing equations for our binary fluid system 
are the Cahn-Hilliard equation for the composition 0, 
and the incompressible Navicr-Stokes equation for the 
velocity Ua in an isothermal fluid of unit density p: 

{dtUa+ UpV fjUa) +V aP - vV'^Ua - (fN alJ- = 0(1) 
dtcj} + Va {<j)U„ - MVa^i) = (2) 

Here, p is pressure (related in LB to density fluctuations, 
which arc small v is the kinematic viscosity; M 

is the ((/)-indepcndent) mobility and p, = B(f) (^cfP — l) — 
kV^0 is the chemical potential. B and k are positive 
constants; the interfacial tension is cr = (8kB/9)1/2 and 
the interfacial width is = (2k/B)1/2 

We solve these equations with an LB algorithm simi- 
lar to that reported in ^2, 15 1. To achieve the necessary 
shear rates, the domain is decomposed blockwise using 
multiple Lees-Edwards sliding periodic boundary condi- 
tions |6[ |l^, chosen so that j^^VyUxdy = ^yi- Al- 
though we neglect thermal fluctuations in our fluid, as 
appropriate for dynamics near a zero-temperature fixed 
point |l7j , a fluctuating local velocity field still arises via 
nonlinear interaction between the order parameter field 
and flow field. To help control errors, we adhered as far 
as possible to previously used parameter values and pro- 
tocols H, However, the sheared 3D case showed sig- 
nificant stability problems compared with either the 3D 
unsheared case or 2D sheared case. To alleviate these, 
we replaced the D3Q15 lattice of with a D3Q19 
model; this removes a "computational mode" responsi- 
ble for some of the instabilities of D3Q15 We also 
use a multiple relaxation time approach p9[ | in place of a 
single relaxation time 12 , further improving stability. 

Most of our 3D production runs were made using sys- 
tem size 512 X 256 x 256, run for t ~ 4 x 10^ time steps. 
Holding other parameters fixed, one finds that if 7 is too 
small, the domain size is large and finite size effects dom- 
inate, whereas if 7 is too large then the domains become 
small on the lattice scale and tend to form a partially 
(or even fully) remixed state with strongly blurred in- 
terfaces. Such remixing could be a real physical effect 
at shear rates so high that the local interfacial structure 
departs strongly from equilibrium, but this happens at 
much lower shear rates in an LB fluid than in a real 
one (where ^0 is much smaller). We therefore reject as 
artefacts all such partially remixed states, as identified 
by a significant reduction in order parameter variance 
(0^). Worst affected were the runs at higher Reynolds 
number (low viscosity) where an adjustment of the in- 
terfacial width from ^0 = 1-13 to ^0 = 1-35 helped to 
maintain acceptable behavior. All simulations reported 
here were done for fully symmetric quenches with param- 
eters summarized in Table 1. As in the unsheared case 
judicious combinations of ^0, c, M and v allow sys- 
tems spanning several decades in L/Lf) and 7T0 to be 
accurately studied by varying Lq ^md Tq alongside 7. 

Fig. |l| shows snapshots of the interfacial structure 
based on the order parameter field for R019 with 7 = 
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R028 


1.41 


0.05 


0.063 


0.055 


36.1 


927 


1.13 


5.0x10" 


-4 


1024 


- 


- 


- 


R029 


0.2 


0.15 


0.047 


0.042 


0.952 


4.54 


1.13 


5.0x10" 


-4 


1024 


- 


- 


- 


R020 


0.025 


2 


0.0047 


0.0042 


0.149 


0.886 


1.13 


5.0x10" 


-4 


1024 


- 


- 


- 


R003 


0.015 


2.0 


0.0047 


0.0042 


0.054 


0.19 


1.13 


7.5x10" 


-4 


512 


511 


72.2 


172 


















5.0x10" 


-4 


512 


828 


116 


352 


R004 


0.01 


2.0 


0.0047 


0.0042 


0.024 


0.0567 


1.13 


7.5x10" 


-4 


512 


356 


68.1 


131 


















5.0x10" 


-4 


512 


491 


106 


192 


R030 


0.00625 


1.25 


0.0047 


0.0042 


0.00930 


0.0138 


1.13 


5.0x10" 


-4 


512 


375 


91.6 


160 


R007 


0.005 


2.0 


0.0047 


0.0042 


0.0059 


0.00709 


1.13 


5.0x10" 


-4 


512 


382 


97.4 


174 


R008 


0.0035 


2.0 


0.0047 


0.0042 


0.0029 


0.00243 


1.35 


5.0x10" 


-4 


512 


370 


101 


177 


R019 


0.0014 


4 


0.0024 


0.0021 


0.000933 


0.000622 


1.35 


5.0x10" 


-4 


512 


234 


71.3 


118 


R032 


0.0005 


5 


0.00094 


0.00083 


0.000301 


0.000181 


1.35 


5.0x10" 


-4 


512 


135 


48.0 


71.2 



Table I: Parameter sets used in 3D simulations and observed NESS length scales. Where a trivial NESS could be identified 
by inspection, no length is recorded. The results of R020 were ambiguous: periods of apparent NESS were contaminated by 
intervals of partial remixing (low {4'^))- 
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Figure 2: Two examples of Lx,y,z in lattice units as a function 
of time in strain units jt for R030 parameters (upper panel) 
and R019 parameters (lower panel). For all parameters where 
a steady state is observed, the length scale as measured by 
the gradient statistic of [^] is largest in the velocity direction 
Lx, followed by the vorticity direction L^, with that in the 
velocity gradient direction Ly the smallest. 

5 X 10"'^ after a steady state had been reached |Q . Fig. |^ 
shows time series for L^^y^z from runs R030 and R019 as 
measured by a standard order parameter gradient statis- 
tic 01 that effectively measures the mean distance be- 
tween interfaces crossing the chosen direction. 

In ||l2| , finite size effects (in the absence of shear) were 
considered quantitatively under control when the correla- 
tion length L was less than 1/4 of the system size A. In 

this criterion was applied to timc-averaged correlation 
lengths Lx^y in the 2D sheared system. However, the ac- 
tual system size dependence of L^^y^z in both 2D § and 
3D (this work) suggests that under shear this criterion is 
unnecessarily strict, at least if the purpose is to eliminate 
the qualitatively artefactual states that arise directly from 
finite size effects. As mentioned previously, these "triv- 
ial" NESS's form obvious laminar stripes extending the 
full size of the simulation box in both x and z directions. 



For such states, L^^z values that are formally much larger 
than the simulation dimensions A^^z are rapidly estab- 
lished. {Lx ^ means that that, for most coordinates 
y, z, one can cycle round the periodic boundary condi- 
tions in X without encountering a single domain wall.) 
To formally eliminate these, a criterion L^^y^z < ^x,y,z 
is applied, which also excludes one apparently nontrivial 
NESS run (Table 1) from the scaling analysis made be- 
low. At the lowest Reynolds numbers investigated, only 
a trivial NESS was found on a 512 x 256 x 256 lattice; 
larger systems, 1024 x 512 x 512, were then simulated for 
these parameters but gave the same structure. This diffi- 
culty in achieving bulk NESS at low Re perhaps suggests 
onset of a new regime; we return to this below. 

Only at the largest {'jTq)^^ values investigated was 
the strict finite size criterion of L^ < A^/^, ap- 
proached. (Note however that earlier studies accepted 
L < A/2 as sufficient, e.g. [^.) Accordingly we ex- 
pect that the quantitative scaling of all our correlation 
length data with shear rate may still be affected by finite 
size corrections. With this caveat, we proceed to per- 
form a scaling analysis based on the protocol of To 
construct our scaling plot, mean values of L^^y^z were ob- 
tained via a bootstrap procedure |^ performed on each 
times series, discarding data for which t < 10^ (to elim- 
inate transients). The results for L^^y^z/Lo are plotted 
against (jTo)^^ in Fig. |3[ Linear least-squares fits to 
these data suggest scaling exponents for L^^y^z of, re- 
spectively, -0.54 ± 0.03, -0.65 ± 0.03, and -6.60 ± 0.04 
at the 95% confidence level. An alternative scaling, us- 
ing the principal axes of the gradient statistic |Q, |j , gives 
exponents for L^,Li,L2 of -0.53 ± 0.04, -0.67 ± 0.03, 
and —0.64 ± 0.06 (data not shown). These results ap- 
pear to rule out Ly ^ 'y~3/4 ^q^^ found in 2D 
However, the range of Re accessible is restricted to 
about 1 decade (260< Re < 2300); as in 2D, one cannot 
rule out that these are effective exponents describing the 
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Figure 3: Reduced length scales Lx,y,z/L[) (black squares, 
red triangles, blue circles respectively) as a function of inverse 
reduced shear rate for the 8 runs in which nontrivial NESS 
was observed. The standard errors in the individual points 
are no larger than the symbols; the dashed lines give the 95% 
confidence limits of the fitted regression. 

crossover region. These Re values are also high enough 
that a multiple length scaling might be needed |2^ . 

The quoted error margins do not, of course, allow for 
systematic error of which there are several sources (even 
discounting finite size effects), each at the likely level of 
several percent ||, ^ . Accordingly these results do not 
rule out a common scaling of all three correlation lengths 
with a single exponent, L^^y^z/Lo ^ {'^To)~^^^ , at least 
in the inertial limit of very large (7ro)^^ where the data 



hints that the three curves may saturate to fixed ratios. 
Conversely, the ever increasing difficulty to achieve NESS 
at small (7To)~^ may point to a quite different behav- 
ior at low Reynolds numbers. Suggestively, Fielding 
has recently performed 2D binary Stokes flow simulations 
finding no evidence of bulk NESS at Re = 0; this could 
mean that inertia plays the role of a singular perturba- 
tion in this problem. Moreover, for a range of 7T0 around 
10"^, NESS is easily achieved in 2D but not 3D: the abil- 
ity to form connections in the vorticity direction might, at 
moderate and low Re, require formation of domains of ex- 
tremely high aspect ratio before a NESS can be reached. 

In conclusion, while open issues remain concerning the 
details of scaling and finite size behavior, our simulations 
present clear evidence for nonequilibrium steady states in 
3D sheared binary fluids. The qualitative character of the 
NESS achieved in these simulations at high Re (low shear 
rate) , which entails a balance between domain stretching 
under flow and coarsening driven by interfacial tension, 
strongly suggests that these results represent true bulk 
behavior. Since the effect of coarsening at flxed 7 is to 
increase Re, indefinite coarsening |^ can seemingly be 
ruled out even at higher shear rates, although a different 
mechanism for achieving NESS may operate there. 
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